In [2]:
import nibabel as nib
import nilearn
from nilearn import image
import numpy as np
import pandas as pd
import os.path as osp
from glob import glob
from nilearn import plotting
In [17]:
c1dir = '/media/Projects/ALFA_VBM_processed/c1/'
mddir = '/home/grg/data/ALFA_DWI/'
jacdir = '/home/grg/spm/Jacobians/'
perdir = '/home/grg/spm/ROIvent'
dm = pd.read_excel('/home/grg/spm/designmatrix3.xls')
vv = pd.read_excel('/home/grg/spm/data/Ventricular volumes.xlsx')
fsvol = pd.read_excel('/home/grg/spm/data/aseg FS ALFA.xlsx')
subjects_fp = '/home/grg/spm/data/subjects.json'
import json
allsubjects = json.load(open(subjects_fp))
ages = json.load(open('/home/grg/spm/data/age.json'))
educyears = dict([(str(int(e)), v) for e, v in dm[['Subj_ID', 'Years of Education']].to_dict(orient='split')['data']])
genders = dict([(str(int(e)), v) for e, v in dm[['Subj_ID', 'Gender(0=female)']].to_dict(orient='split')['data']])
tivs = dict([(str(int(str(e)[:5])), v) for e, v in vv[['subject', 'Total Intracranial Volume']].to_dict(orient='split')['data']])
c1fp = '/media/Projects/ALFA_VBM_processed/c1/c1t1_10013.nii'
perfp = osp.join(perdir, '10013_latvent_dilated.nii.gz')
fsvol = fsvol.loc[fsvol['StructName'].isin(['Left-Lateral-Ventricle', 'Right-Lateral-Ventricle', 'Left-Inf-Lat-Vent', 'Right-Inf-Lat-Vent'])][['subject', 'Volume_mm3', 'StructName']]
dm[dm['Subj_ID'] == 10070]
Out[17]:
In [ ]:
In [ ]:
In [18]:
data = []
subjects = allsubjects #[0:100]
groups = ['Apoe2-2', 'Apoe2-3', 'Apoe2-4', 'Apoe3-3', 'Apoe3-4', 'Apoe4-4']
agegroups = ['Apoe2-3','Apoe2-4','Apoe3-3','Apoe3-4','Apoe4-4','age23',
'age24','age33','age34','age44','agesq23','agesq24','agesq33','agesq34','agesq44']
for j, s in enumerate(subjects):
print j, s
#try:
mdfp = glob(osp.join(mddir, '%s*'%s, 'DWI', '%s*_MD_t1space.nii.gz'%s))[0]
jacfp = glob(osp.join(jacdir, '*%s*nii'%s))[0]
perfp = glob(osp.join(perdir, '%s*.nii.gz'%s))[0]
md = np.array(nib.load(mdfp).dataobj, copy=True)
peri = np.array(nib.load(perfp).dataobj, copy=True)
mdl = np.mean(md[peri==1])
mdr = np.mean(md[peri==2])
vvol = dm[dm['Subj_ID'] == s]['ventricles_FS'].tolist()[0]
age = ages[str(s)]
agesq = age * age
tiv = tivs[str(s)]
sex = genders[str(s)]
ey = educyears[str(s)]
#for group, g in enumerate(groups):
# if s in dm[dm[g] == 1]['Subj_ID'].tolist():
# break
#apo = [dm[dm['Subj_ID'] == s][x].tolist()[0] for x in agegroups]
data.append([mdl, mdr, vvol, age, agesq, sex, ey]) #, tiv, group])
#data[-1].extend(apo)
#except:
# print 'skipping', s
In [19]:
columns = ['mdl', 'mdr', 'vvol', 'age', 'agesq', 'gender', 'educyears']
#columns.extend([e.replace('-','') for e in agegroups])
data
Out[19]:
In [440]:
print apo
df = pd.DataFrame(data, columns=columns)
df
Out[440]:
In [416]:
from pandas.tools.plotting import scatter_matrix
%matplotlib inline
import matplotlib.pyplot as plt
axes = scatter_matrix(df, figsize=(35,35))
corr = df.corr().as_matrix()
for i, j in zip(*plt.np.triu_indices_from(axes, k=1)):
axes[i, j].annotate("%.3f" %corr[i,j], (0.8, 0.8), xycoords='axes fraction', ha='center', va='center')
plt.show()
In [16]:
from scipy import stats
import patsy
import statsmodels.api as sm
from statsmodels.formula.api import ols
print columns
f= 'md_ofcl ~ C(apoe) + age:C(apoe) + agesq:C(apoe) + C(gender) + educyears'
lm = ols( f, #Apoe23 + Apoe24 + Apoe33 + Apoe34 + Apoe44 + age23 + age24 + age33 + age34 + age44 + \
#agesq23 + agesq24 + agesq33 + agesq34 + agesq44 + gender + educyears',
data=df).fit()
y,X = patsy.dmatrices(f, df, return_type='dataframe')
table = sm.stats.anova_lm(lm)
r = np.zeros_like(lm.params)
r[1] = 3
r[2]= -2
r[3] = 3
r[4] = -2
r[5] = -2
print lm.params
hypothese = 'C(apoe)[T.5] - (C(apoe)[T.4] + C(apoe)[T.3] + C(apoe)[T.2])'
lm.t_test(hypothese)
In [403]:
box = [df[df['apoe'] == x]['md_ofcl'].tolist() for x in xrange(6)]
plt.boxplot(box)
plt.show(box)
In [391]:
from nistats.thresholding import map_threshold
im = np.array(nib.load('/tmp/ofc_mask2.nii.gz').dataobj, copy=True)
thresholded_map1, threshold1 = map_threshold('/tmp/output_FSvent/MD/analysis/estimatecontrasts/spmT_0006.nii', threshold=0.001, cluster_threshold=50)
im2 = np.array(thresholded_map1.dataobj, copy=True)
im[im2==0] =0
from nilearn import plotting
plotting.plot_glass_brain(image.new_img_like('/tmp/ofc_mask.nii.gz', im))
image.new_img_like('/tmp/ofc_mask.nii.gz', im).to_filename('/tmp/ofc_mask2_shrinked.nii.gz')
In [19]:
from skimage import morphology as morpho
from skimage.morphology import ball
import numpy as np
import nibabel as nib
from nilearn import image
from nilearn import plotting
import os.path as osp
%matplotlib inline
%run /home/grg/git/alfa/nilearn-helper.py
In [52]:
vent_fp = '/home/grg/data/ALFA_DWI/10013/T1/10013_latvent.nii'
#t1_fp = '/home/grg/data/ALFA_DWI/10013/T1/10013_mabonlm_nobias.nii'
a1 = np.array(nib.load(vent_fp).dataobj)
a2 = morpho.dilation(a1, ball(6))
mit = a2.shape[0] /2
print mit
a2[a2<1] = 0
a2[mit:,:,:][a2[mit:,:,:] > 0.5] = 2
a2[a1 != 0] = 0
image.new_img_like(vent_fp, a2).to_filename('/tmp/test.nii.gz')
In [30]:
img = plotting.plot_roi(image.new_img_like(vent_fp, a), bg_img=t1_fp, display_mode='x', cut_coords=range(-10,10,2))
#dd_overlay(image.new_img_like(vent_fp, a))